#include <stdlib.h>
#include <stdio.h>
#include "nbody.h"
#include <math.h>
#include <string.h>
#include <omp.h>

#define MAX 255

int numCuerpos, numPasos, anchuraTablero, alturaTablero,numProcesos;
double g, intervalo, minDist;
double *posicion[2];
double *velocidad[2];
double *fuerza[2];
double *masa;
int *indices[4];


int main(int argc, char **argv){
    int i;
    leerCabecera(argv[1]);
    numProcesos=atoi(argv[3]);
    omp_set_num_threads(numProcesos);
    calcularIndices();
    for(i=0;i<numPasos;i++){
        calcularNbody();
    }
    
    printf("Guardando los datos en el archivo: %s\n", argv[2]);
    escribirDatos(argv[2]);
    
    return 1;
}

int leerCabecera(char *fichero_entrada){
    FILE *file;
    int cont;
    
    //Apertura del fichero y lectura de la cabecera
    
    file = fopen(fichero_entrada, "r");
    
    fscanf(file, "%d %lf %lf %d %lf %d %d", &numCuerpos, &g, &intervalo, &numPasos, &minDist, &anchuraTablero, &alturaTablero);
    
    inicializarDatos();
    
    printf("numcuerpos: %d, g: %lf, intervalo: %lf, numPasos: %d,mindist: %lf,anchura: %d,altura: %d\n", numCuerpos, g, intervalo, numPasos, minDist, anchuraTablero, alturaTablero);
    
    for(cont=0; cont<numCuerpos; cont++){
        fscanf(file, "%lf %lf %lf %lf %lf \n", &posicion[0][cont], &posicion[1][cont], &velocidad[0][cont], &velocidad[1][cont], &masa[cont]);
    }
    return 1;
}

int escribirDatos(char *archivo){
    FILE *file;
    int cont;
    
    //Apertura del fichero y escritura de la cabecera
    file = fopen(archivo, "wr");
    fprintf(file, "%d %lf %lf %d %lf %d %d\n", numCuerpos, g, intervalo, numPasos, minDist, anchuraTablero, alturaTablero);
    
    //Escritura de los datos
    for(cont=0; cont<numCuerpos; cont++){
        fprintf(file, "%lf %lf %lf %lf %lf\n", posicion[0][cont], posicion[1][cont], velocidad[0][cont], velocidad[1][cont], masa[cont]);
    }
    fclose(file);
    return 1;
}

void inicializarDatos(){
    posicion[0] = (double *)malloc(numCuerpos* sizeof (double));
    posicion[1] = (double *)malloc(numCuerpos* sizeof (double));
    fuerza[0] = (double *)malloc(numCuerpos* sizeof (double));
    fuerza[1] = (double *)malloc(numCuerpos* sizeof (double));
    velocidad[0] = (double *)malloc(numCuerpos* sizeof (double));
    velocidad[1] = (double *)malloc(numCuerpos* sizeof (double));
    masa = (double *)malloc(numCuerpos* sizeof (double));
    indices[0]=(int *)malloc(numProcesos*sizeof(int));
    indices[1]=(int *)malloc(numProcesos*sizeof(int));
    indices[2]=(int *)malloc(numProcesos*sizeof(int));
    indices[3]=(int *)malloc(numProcesos*sizeof(int));
}

void calcularNbody(){
    double fax, fay, dist, f, alpha, aux, *ax, *ay;
    int i, j, fin=0;
    
    
    #pragma omp parallel for
    for (i = 0; i < numCuerpos; i++) {
        fuerza[0][i] = fuerza[1][i] = 0;
    }
    
    #pragma omp parallel firstprivate(fin) private(i,j,dist,f,aux,alpha,fax,fay) 
    {
    for (i = indices[0][omp_get_thread_num()]; (i <= indices[2][omp_get_thread_num()]&& !fin); i++) {
        
        if (i== indices[0][omp_get_thread_num()])
            j=  indices[1][omp_get_thread_num()];
        else
            j=  i+1;
        
        for (; j<numCuerpos && !fin; j++) {
            dist = sqrt((posicion[0][j] - posicion[0][i]) * (posicion[0][j] - posicion[0][i]) + (posicion[1][j] - posicion[1][i]) * (posicion[1][j] - posicion[1][i]));
            
            if (dist > minDist ) {
                
                f = g * masa[i] * masa[j] / (dist * dist);
                
                if (f > MAX_F)
                    f = MAX_F;
                
                aux  = (posicion[1][j]-posicion[1][i]) / (posicion[0][j]-posicion[0][i]);
                
                alpha = atan(aux);
                
                fax   = f * mycos(alpha);
                fay   = f * mysin(alpha);
                
#pragma omp critical
                {
                    fuerza[0][i] += fax;
                    fuerza[1][i] += fay;
                    fuerza[0][j] -= fax;
                    fuerza[1][j] -= fay;
                }
                
            }
            if((j+1)==indices[3][omp_get_thread_num()]&&i == indices[2][omp_get_thread_num()])
                fin=1;
        }
    }
    }

    
    ax = (double *)malloc(sizeof(double) * numCuerpos);
    ay = (double *)malloc(sizeof(double) * numCuerpos);
    
     #pragma omp parallel for
    for (i = 0; i < numCuerpos; i++) {
        
        //calculo de las aceleraciones
        ax[i] = fuerza[0][i] / masa[i];
        ay[i] = fuerza[1][i] / masa[i];
        //calculo de las velocidades
        velocidad[0][i] = velocidad[0][i] + ax[i] * intervalo;
        velocidad[1][i] = velocidad[1][i] + ay[i] * intervalo;
        //calculo de las posiciones
        posicion[0][i] = posicion[0][i] + velocidad[0][i] * intervalo;
        posicion[1][i] = posicion[1][i] + velocidad[1][i] * intervalo;
        // printf("[%d] px: %lf vx: %lf intervalo: %lf\n",i,posicion[0][i],velocidad[0][i],intervalo);
        if (posicion[0][i] <= 0) {
            posicion[0][i] = 2;
            velocidad[0][i] = - velocidad[0][i];
        }
        if (posicion[0][i] >= anchuraTablero ) {
            posicion[0][i] = anchuraTablero - 2;
            velocidad[0][i] = - velocidad[0][i];
        }
        if (posicion[1][i] <= 0) {
            posicion[1][i] = 2;
            velocidad[1][i] = - velocidad[1][i];
        }
        if (posicion[1][i] >= alturaTablero) {
            posicion[1][i] = alturaTablero - 2;
            velocidad[1][i] = - velocidad[1][i];
        }
        //printf("body %d: x=%d, y=%d vx=%lf vy=%lf ax=%lf ay=%lf fx=%lf  fy=%lf\n", i, x[i], y[i],vx[i],vy[i],ax[i],ay[i],fx[i],fy[i]);
    }
    
}


void calcularIndices(){
    
    int numIteracionesTotales=0, numIteProceso, i, j, aux=0, procesoActual=1;
    
    for(i=1; i<numCuerpos; i++)
        numIteracionesTotales+=i;
    
    numIteProceso=numIteracionesTotales/numProcesos;
    
    indices[0][0]=0;
    indices[1][0]=1;
    
    for (i=0;i<numCuerpos&&procesoActual<(numProcesos);i++) {
        for(j=i+1;j<numCuerpos&&procesoActual<(numProcesos);j++) {
            aux++;
            if(aux==numIteProceso) {
                
                if((j+1)<numCuerpos) {
                    indices[2][procesoActual-1]=i;
                    indices[3][procesoActual-1]=j+1;
                    indices[0][procesoActual]=i;
                    indices[1][procesoActual++]=j+1;
                }
                else {
                    indices[2][procesoActual-1]=i+1;
                    indices[3][procesoActual-1]=i+2;
                    indices[0][procesoActual]=i+1;
                    indices[1][procesoActual++]=i+2;
                }
                aux=0;
            }
        }
    }
    
    if((j)<numCuerpos) {
        indices[0][numProcesos-1]=i-1;
        indices[1][numProcesos-1]=j;
    }
    else {
        indices[0][numProcesos-1]=i-1;
        indices[1][numProcesos-1]=i+1;
    }
    indices[2][numProcesos-1]=numCuerpos;
    indices[3][numProcesos-1]=numCuerpos;
}
double factorial(double x) {
    double fac = 1;
    int i;
    for (i=2; i<=x; i++) fac *= i;
    return fac;
}

double mypow(double x, double exp) {
    double result = 1;
    int i;
    
    result = 1;
    for (i=0; i<exp; i++) result *= x;
    return result;
}

double mysin(double alpha) {
    double result = 0;
    int i;
    
    for (i=0; i < MAX_ACCURACY;i++) {
        result += (double)(mypow(-1, i) / factorial((double)(2*i+1))) * (double)mypow(alpha, 2*i+1);
    }
    return result;
}

double mycos(double alpha) {
    double result = 0;
    int i;
    
    for (i=0; i < MAX_ACCURACY;i++) {
        result += (double)(mypow(-1, i) / factorial((double)(2*i))) * (double)mypow(alpha, 2*i);
    }
    return result;
}

